Phase behaviour and particle-size cutoff effects in polydisperse fluids 
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We report a joint simulation and theoretical study of the liquid- vapor phase behaviour of a fluid 
in which polydispersity in the particle size couples to the strength of the interparticle interactions. 
Attention is focussed on the case in which the particles diameters are distributed according to a 
fixed Schulz form with degree of polydispersity 5 = 14%. The coexistence properties of this model 
are studied using grand canonical ensemble Monte Carlo simulations and moment free energy calcu- 
lations. We obtain the cloud and shadow curves as well as the daughter phase density distributions 
and fractional volumes along selected isothermal dilution lines. In contrast to the case of size- 
mdependent interaction strengths (N.B. Wilding, M. Fasolo and P. Sollich, J. Chem. Phys. 121, 
6887 (2004)), the cloud and shadow curves are found to be well separated, with the critical point 
lying significantly below the cloud curve maximum. For densities below the critical value, we ob- 
serve that the phase behaviour is highly sensitive to the choice of upper cutoff on the particle size 
distribution. We elucidate the origins of this effect in terms of extremely pronounced fractionation 
effects and discuss the likely appearance of new phases in the limit of very large values of the cutoff. 
PACS numbers: 64.70Fx, 68.35.Rh 



I. INTRODUCTION AND BACKGROUND 



When the constituent particles of a many body sys- 
tem exhibit variation in terms of some attribute such as 
their size, shape or charge, then that system is termed 
"polydisperse" . Occurring, as it does, in soft matter sys- 
tems as diverse as colloids, polymers and liquid crystals, 
the phenomenon of polydispersity is both widespread and 
basic. However, owing to the complexity that polydisper- 
sity bestows on a system, gaining an understanding of its 
physical consequences represents a considerable challenge 
to theory, simulation and experiment alike. Meeting this 
challenge is not solely a matter of fundamental concern, 
but also one of practical importance: the presence of 
polydispersity in substances of commercial and indus- 
trial importance (such as paints, fuels and foodstuffs) 
is known to profoundly affect their thermodynamic and 
processing properties in ways which are as __^et neither 
well characterized nor well understood 



ire as ye 
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Of the many fundamental questions that polydisper- 
sity elicits, one of the most central concerns its influ- 
ence on bulk phase separation. The phase behavior of 
polydisperse fluids is known to be considerably richer in 
both variety and character than that of their monodis- 
perse counterparts |3lPL This richness has its source in 
fractionation effects [ill B H E IH IH US : the dis- 
tribution of the polydisperse attribute differs from one 
coexisting phase to another. In order to quantify this 
effect, it is expedient Q to regard the polydisperse at- 
tribute as a continuous variable a, and to define a density 
distribution p{cr), where p(a)da is the number density of 
particles in the range a . . . a + da. Fractionation then oc- 
curs when, at two phase coexistence, a "parent" density 
distribution p^^'{a) splits into two distinct "daughter" 
phase distributions ^'""^'((7) and p*^^-'(ct). Particle conser- 
vation implies that the daughter distributions are related 



to the parent via a simple volumetric average: 



il-OP^'\T)+^p^^Ha) 



^P^"H-), 



(1) 



where 1 — ^ and ^ are the fractional volumes of the two 
phases. 

Experimentally, the situation typically considered (for 
example in studies of colloidal dispersions) is the nature 
of the phase behaviour along a so-called dilution line. 
Here, one takes a sample of some prescribed polydisper- 
sity and observes its coexistence properties as it is pro- 
gressively diluted by adding solvent, the temperature be- 
ing held constant. Under such circumstances the parent 
distribution maintains a fixed shape, and only its overall 
scale varies according to the degree of dilution. Accord- 
ingly, one can write 



pW(a)=n(«)/W(a), 



(2) 



where f^^' (a) is a fixed shape function which serves to de- 
fine the form of the polydispersity, while n^^' is the over- 
all (parent) number density of the system, whose value 
parameterizes the location of the system on the dilution 
line. 

In order to construct the full phase diagram for such a 
system one must obtain the dilution line phase behaviour 
for the range of temperatures of interest. However, in 
contrast to a monodisperse system where the limits of 
phase stability and the densities of the coexisting phases 
are completely specified by the coexistence binodal in 
the n^^' — T plane, fractionation implies that polydis- 
persity splits the binodal into cloud and shadow curves 
,_5j , as shown schematically in fig. ^ These curves mark, 
respectively, the density of the onset of phase separation 
and the density of the incipient (shadow) daughter phase. 
For parent densities lying wholly within the coexistence 
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FIG. 1: A schematic fluid-fluid phase diagram for a polydis- 
perse fluid, showing temperature T against density; the cloud 
curve gives the parent density rr ' where phase separation 
first occurs while the shadow curve records the density of the 
incipient coexisting phases. Note, however, that the shadow 
phase density distribution lies off the dilution line. 



region, two daughter phases form, the properties of which 
vary non-trivially with n^'^\ Hence a full specification of 
the coexistence properties of a polydisperse system re- 
quires not only a determination of the locus of the cloud 
curve, but also the dependence of ^, p'^-'^^(cr) and p'^'^\a) 
on ny'' and T . 

In the present work, we describe a joint simulation and 
theoretical study of liquid-vapor phase separation for a 
model fluid in which polydispersity affects not only the 
length scale but also the strength of the interaction po- 
tential. Previously we have considered the restricted case 
of size-independent interaction strengths, for which the 
critical point is found to lie very close to the maximum of 
the cloud curve Tifl. By contrast, for the present model it 
is located substantially below the maximum (cf. fig. Q), 
as is observed in many experiments on complex fluids 
(see e.g. [ia|)- One implication of this is that phase co- 
existence occurs at and above the critical temperature, 
provided the overall parent density is less than its critical 
value, and we investigate this feature here. Another dif- 
ference to the case of size-independent interactions is an 
acute sensitivity of the coexistence properties to the pres- 
ence of rare large particles. We study this phenomenon 
by varying a particle size cutoff parameter. 

As regards previous simulation studies of polydisperse 
phase equilibria, other authors have (to date) considered 
exclusively the case of variable polydispersity, in which 
the shape of the parent distribution chan ges a ccording to 
the temperature and overall density [l^ [13, llSi J^: ^^ 
l2l| . As such, these studies do not reflect the most com- 
monly encountered experimental situation in which the 
polydispersity of a system is fixed by the process of its 
chemical synthesis. With the recent advent of bespoke 
Monte Carlo (MC) simulation techniques [i3j|22|j how- 
ever, studies of fixed polydispersity within the (appropri- 
ate) grand or semi-grand canonical ensemble framework 
are now tractable. In the present work we apply these 



methods to study the phase separation of a fluid whose 
polydispersity assumes a prescribed functional form. 

A number of analytical studies of polydisperse phase 
equilibria have also been reported in the literature. These 
typically seek to calculate the system free energy as a 
function of a set of density variables. Unfortunately, this 
task is complicated by the fact that the requisite free en- 
ergy f[p'''^\o)] is a functional of p^^\a), and therefore 
occupies an infinite dimensional space. For sufficiently 
simple model free energies which generalize the van der 
Waals (vdW) approximation, a direct attack on the so- 
lution of the phase equilibrium conditions is nevertheless 
sometimes possible, see e.g. 0i|23- The reason for this is 
that such models are normally "truncatable" [2J| so that 
the phase equilibrium conditions can be reduced to non- 
linear equations for a finite number of variables. This 
approach has been applied to the study of phase sep- 
aration in fluids exhibiting separate size and interaction 
strength polydispersity, yielding predictions for the cloud 
and shadow curves and critical parameters as a function 
of polydispersity |25|. An alternative approach which 
more systematically exploits the advantages of truncat- 
able models is the moment free energy (MFE) method 
1^ |24| . This approximates the full free energy appropri- 
ately in terms of a "moment free energy" which depends 
on a small number of density variables, thereby permit- 
ting the efficient use of standard tangent plane construc- 
tion to locate phase boundaries. The MFE method de- 
livers (for the given model free energy) exact results for 
the location of critical points and the cloud and shadow 
curves. With appropriate refinements, it can also be used 
to obtain exact numerical solutions to the phase coexis- 
tence conditions when two or more phases coexist in fi- 
nite amounts HHHIillill. The MFE method has been 
applied to the study of phase behaviour in systems rang- 
ing from polydisperse hard rods to the freezing of hard 
spheres 0, |23, IM, |3ll |32| , and we shall deploy it again 
here to study the phase behaviour of the present model, 
based on a refined van der Waals-type approximation to 
the excess free energy. 

Our paper is organized as follows [S^l- In sec. ^ we in- 
troduce our model, a polydisperse Lennard-Jones (LJ) 
fluid that we study via grand canonical Monte Carlo 
(MC) simulation and the MFE method. Additionally 
we sketch the MC simulation methodology and provide 
some brief background to the MFE calculations. Sec. IIIII 
details our estimates for the cloud and shadow curves, 
as well as the behaviour of the fractional volumes of the 
daughter phases and their density distributions along se- 
lected isothermal dilution lines in a range including the 
critical temperature. We then turn in sec. lIVI to a consid- 
eration of the effects of the choice of the upper cutoff on 
the particle size distribution. Our results demonstrate 
that coexistence properties can be acutely sensitive to 
this choice. We elucidate the origin of this sensitivity in 
terms of the interplay between the rate of decay of the 
parent size distribution at large a^ and the cr-dependence 
of the excess chemical potential. Likely implications for 



the phase behaviour at very large cutoffs are also dis- 
cussed. Finally, Sec. [V] details out conclusions and high- 
lights some outstanding questions worthy of further in- 
vestigation. 



II. MODEL AND METHODOLOGIES 

A. Simulation model 

We consider a polydisperse Lennard-Jones fluid whose 
interparticle potential takes the form: 



(cTyMj)^ - {(^ij/nj) 



(3) 



with interaction strength e^j = UiUj and interaction ra- 
dius CTij ~ {(7i + o-j)/2; rij = |ri — Tjl denotes the sep- 
aration between the particles. A truncation was applied 
to the potential for r^ > 2.5cry and no tail corrections 
were used. 

With regard to the form of this potential, a couple 
of remarks are appropriate. Firstly, although it differs 
from that studied in Ref. |l3l | only in the cr-dependence 
of the interaction strength parameter e^j , we shall show 
that this difference is crucial in determining many of the 
qualitative aspects of the phase behavior of the model. 
Secondly, the form of Eq. © might, in one sense, be 
regarded as somewhat artificial because for vanishingly 
small particles the interaction strength approaches zero. 
An arguably more realistic form would be a repulsive 
hard core plus an attractive well of variable strength. 
However, it transpires that for the form of parent dis- 
tribution studied in the present work this drawback is 
of minor significance since the physics of the system is 
dominated by the largest particles. 



B. Parent distribution 

Throughout this work we consider the case in which 
the particle diameters a are drawn from a (parental) dis- 
tribution of the Schulz form ^Jl : 



f^'\o) 



1 



1 



z+l 



a exp 



z + \ 



(4) 



with a mean diameter a which sets our unit length scale. 
We have elected to study the case z = 50, correspond- 
ing to a moderate degree of polydispersity: the standard 
deviation of particle sizes is (5 = I/a/z + 1 w 14% of the 
mean. The form of the distribution is shown in Fig. |21 
Although our motivation for employing the Schulz dis- 
tribution is primarily ad-hoc, we note that it has been 
found to fairly accuratel y d escribe the polydispersity of 
some polymeric systems p35| . 

In both the simulations and the MFE calculations de- 
scribed below, the distribution /'■'^•'((t) was limited to 
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FIG. 2: Schulz distribution for the case z = 50 (cf. Eq. @) 



within the range 0.5 < a < ac- The upper cutoff CTc 
serves to prevent the appearance of arbitrarily large par- 
ticles in the simulation, but would also be expected in 
experiment because in the chemical synthesis of colloid 
particles, time or solute limits restrict the maximum par- 
ticle size [33. Various choices of CTc have been consid- 
ered, and these are discussed in relation to the results of 
Sees, rrm and ITVI 



C. Simulation strategy 

The phase behaviour of the polydisperse LJ model 
Eq. O was studied via MC simulations within the grand 
canonical ensemble (GCE). The MC algorithm invokes 
four types of operation: particle displacements, dele- 
tions, insertions, and resizing. The particle diameter a is 
treated as a continuous variable in the range < cr < ctc , 
with (7c the prescribed cutoff. However, distributions de- 
fined on a such as the instantaneous density p{a) , and the 
chemical potential /x(ct), are represented as histograms 
defined by discretising the permitted range of a into 120 
bins. Most of the simulation results presented below 
pertain to a periodic cubic system of linear dimension 
L = 15ct, although near the critical point, where finite- 
size effects are important, we determined the phase be- 
haviour using system sizes ranging up to L = 30a. For 
further details concerning the simulation algorithm, as 
well as the structure, storage and acquisition of data, we 
refer the interested reader to ref. |33 . 

The principal observable of interest is the fluctuat- 
ing form of the instantaneous density distribution p{a). 
From this we derive the distribution p{n) of the overall 
number density n = J dap{<7), and that of the volume 
fraction 7/ = (tt/G) J daa^p{a). The existence of phase 
coexistence at given chemical potentials is signalled by 
the presence of two distinct peaks in the probability dis- 
tribution p{n). In order to obtain estimates of dilution 
line coexistence properties at some prescribed tempera- 
ture, we employ an approach recently proposed by our- 
selves in Ref. JS^- For a given choice of n^"' , the method 
entails tuning the chemical potential distribution p{a) 



together with a parameter ^, such as to simultaneously 
satisfy both a generalized lever rule and an equal peak 
weight criterion j^ for p{n) : 



,(0) 



r = 1 (5b) 



Here the daughter density distributions p'^^\a) and 
p^'^'{a) are assigned by averaging only over configura- 
tions belonging to either peak of pi^n) . The quantity r is 
the peak weight ratio: 






(6) 



with n* a convenient threshold density intermediate be- 
tween vapor and liquid densities, which we take to be the 
location of the minimmn in p{n). The tuning of fJ.{a) and 
^ necessary to simultaneously satisfy Eqs. H5a|) and (|5b|l 
can be efficiently achieved by a combination of histogram 
extrapolation techniques |43| , and a non-equilibrium po- 
tential refinement (NEPR) procedure as described else- 
where |22j. 

The value of £, resulting from the application of the 
above procedure is the desired fractional volume of the 
liquid phase at the nominated value of n^^\ Cloud 
points are determined as the value of n^'^' at which ^ 
first reaches zero (vapor branch) or unity (liquid branch) , 
while shadow points are given by the density of the co- 
existing incipient daughter phase, which may be simply 
read off from the appropriate peak density in the cloud 
point form of p{n). It should be pointed out that the 
finite-size corrections to estimates of coexistence proper- 
ties obtained using the equal peak weight criterion for 
p{n) are exponentially small in the system size |33| . 

In order to obtain the phase behaviour of our model 
system, we scanned the dilution line for a selection of 
fixed temperatures. We started by setting T — Tc, the 
critical temperature, and tracked the locus of the dilution 
line in a stepwise fashion. This tracking procedure must 
be bootstrapped with knowledge of the form of /i(cr) at 
some initial point on the dilution line. A suitable esti- 
mate was obtained, for a point near the critical density, 
by means of the NEPR procedure |22|, in conjunction 
with the equal peak weight criterion for p{n), discussed 
above. Simulation data accumulated at this near-critical 
state point was then extrapolated to a lower, but nearby 
density n'"' by means of histogram reweighting, thus pro- 
viding an estimate of the corresponding form of /i((t). 
The latter was employed in a new simulation, the results 
of which were similarly extrapolated to a still lower value 
of n^^\ Iterating this procedure thus enabled the sys- 
tematic tracking of the whole dilution line. Histogram 
extrapolation further permitted a determination of di- 
lution line properties at adjacent temperatures, thereby 
facilitating a systematic determination of the phase be- 
haviour in the n^'^' — T plane. Implementation of mul- 
ticanonical preweighting techniques |4lj at each coexis- 



tence state point ensured adequate sampling of the co- 
existing phases in cases where they are separated by a 
large interfacial free energy barrier (see refs. [Ij, |43l fo'^ 
a fuller account of this latter procedure) . 



D. Moment free energy method 

We complement the simulations with theoretical phase 
behaviour calculations, following closely our study of 
the purely size-polydisperse case |l^. To find a suit- 
able expression for the excess free energy density f^^, 
we approximate the repulsive part of our LJ interaction 
as completely hard. The resulting contribution to f^"^ 
is represented accurately by the BMCSL approximation 
/^MCSL Ullillii: 



■^/3/bmcsl 



|-„)mi- 



A3J 



pI 



3piP2 ^ 

1-P3 P3(l-P3)^ 



This is a function of the moments up to third order 
of the density distribution, pi = (jr/G) Jda a^p{a) {i — 
0, . . . , 3); note that po = (7r/6)n is proportional to the 
overall number density, while pa = 77 is the volume frac- 
tion. We have set fee = 1 and defined (3 — 1/T as usual 
so that Pf^^ has dimension of density. 

We treat the remaining attractive part of the LJ inter- 
action potential in the simplest possible way, by adding 
a quadratic van der Waals term to the excess free energy. 
Using the fact that the interaction volume of two parti- 
cles of diameters a and a' scales as {a + a')^, while the 
interaction strength e^ is proportional to ctct', this can 
be written as 



7/3/ 



CX 

vdW 



^ (^) ' J da da' p{a)p{a'){aa'){a + a'f 
-{piPi + ip2P?.) (7) 



where t is an appropriate dimensionless temperature. 
Given the approximate character of the overall excess 
free energy p"" = /^^cgL + /vdw i* ^o^ld not make 
sense to try to scale t precisely to the temperature in the 
simulations. Instead, we will be content to study whether 
our f''^ can reproduce the qualitative trends observed in 
the simulations. While still somewhat crude, it should 
be better suited to this task than previous versions |25l | 
because it incorporates polydispersity not only into the 
attractive contribution, but also into the hard core refer- 
ence system. 

As pointed out in the introduction, to predict phase 
equilibria for polydisperse systems is generally still a 
computationally demanding task even if an explicit ex- 
pression for the excess free energy is available 5]. The 
above excess free energy has the simplifying feature that 
it only depends on the finite set of moments po, . . . ,P4, 
i.e. it is truncatable [2J|. This allows us to exploit the 
MFE method to obtain accurate numerical predictions 
for the phase behaviour. We refer the interested reader 



to our analogous study of the purely size-polydisperse 
case for a fuller description of the methodology [13 • 



III. PHASE BEHAVIOUR: GENERAL ASPECTS 

We have employed the simulation and MFE methods 
outlined in the previous section to determine the cloud 
and shadow curves for our system. Various choices of the 
particle size cutoff have been considered, although the 
majority of our simulation results pertain to the cases 
(Tc = 1.4 and CTc — 1-6. Referring to Fig.|21it is clear that 
both values are quite far out into the tail of the parent 
size distribution, and naively one would therefore expect 
only very minor differences. Our results, expressed in 
the n-T representation, are shown in Fig. Oa) (for sim- 
ulations of the LJ model) and Fig. ^h) (for MFE pre- 
dictions). In both cases we observe a strong separation 
of cloud and shadow curves; as a consequence the criti- 
cal point lies well below the cloud curve maximum. We 
also note a pronounced difference in the cloud curves for 
the two cutoffs, contrary to naive expectation; detailed 
discussion of this effect is deferred to Sec. IIVI 

Finite-size scaling methods |43 (specifically the match- 
ing oip{n) to its known universal critical point form) were 
utilized to provide accurate estimates of the critical pa- 
rameters of the LJ model. These yielded ric — 0.326(3), 



Tc = 1.375(2) (for ac = 1.4); n^ 



(0) 



0.329(3), Tc 



1.384(2) (for ac = 1.6), which are to be compared with 
those of the monodisperse LJ fluid |4a|: ric = 0.3196(4), 
Tc = 1.1876(3). Thus polydispersity of the form consid- 
ered here acts to raise the critical temperature substan- 
tially, although the associated increase in the critical den- 
sity is much more modest. We note that a polydispersity- 
induced increase in Tc has also been observed for the case 
of size-independent interaction strengths [13 ■ There, 
however, the magnitude of the shift in Tc was only around 
5% of that seen in the present model, for the same parent 
distribution. 

The MFE results of Fig. Olb) show qualitatively very 
similar behaviour to the simulation results: the high 
density part of the shadow curve has an unusual posi- 
tive slope which becomes more pronounced as the cutoff 
increases; at the same time the maxima of cloud and 
shadow curves shift to higher temperatures. Quantita- 
tively the cutoff effects are somewhat weaker than in the 
simulations, which is why we chose to display results for 
tJc = 1-6 and 1.9 rather than 1.4 and 1.6 as for the sim- 
ulations. 

One difference between simulations and theory is the 
nature of the cloud and shadow curves in the vicinity 
of the critical point. In the simulation results "dips" 
are observable in this region, in contrast to the MFE re- 
sults which show no such effect. The precise origin of the 
dips is not clear to us at present. One possibility is that 
they are a finite-size effect caused by the breakdown near 
the critical point of our procedure for determining cloud 
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FIG. 3: Cloud and shadow curves in the n — T plane, (a) MC 
simulation results, (b) MFE predictions. Lines are a guide 
to the eye. 



points. Alternatively they may be indicative of genuine 
critical point singularities in the cloud and shadow curves 
which are not picked up by the theory because our ap- 
proximate free energy expression is of a mean-field type. 

In Fig. 01 we show the forms of the cloud and shadow 
curves in the rj-T representation, i.e. plotting on the x- 
axis the volume fraction i] rather than the density n of 
the cloud and shadow phases. Since the cloud phase al- 
ways has the fixed parental size distribution, the cloud 
curve is simply rescaled by this change of representation; 
the same is not true of the shadow curve, however, since 
the size distribution in the shadow varies throughout the 
phase diagram. As a consequence, we observe that in the 
volume fraction plot the cloud and shadow curves sepa- 
rate even further than before, forming a rather symmet- 
rical "butterfly" shape. Again there is good qualitative 
agreement with the MFE predictions. 

Turning now to the coexistence properties inside the 
coexistence region, we plot, in Figs. El and the mea- 
sured n^^' dependence of the number densities and vol- 
ume fractions of the coexisting vapor and liquid phases on 
the critical isotherm. In the MC simulations, these were 
obtained simply by reading off the peak positions of p{n) 
and p{ri) for the appropriate values of n^^' and T. Close 
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FIG. 4: Cloud and shadow curves in the rj — T plane, (a) 
MC results, (b) MFE predictions. Lines are a guide to the 
eye. 



to the critical point, one expects large finite-size effects to 
occur due to the divergent correlation length. These are 
indeed apparent in our measurements of the peak densi- 
ties of p{n) and p{ri) for a number of system sizes (see 
insets of Figs.|S{a) andEJ^a)). Specifically, since the form 
oip{n) is doubly peaked for finite-sized systems, even at 
the critical point [43, one expects that the peak den- 
sities approach the critical density as a universal power 
law in the system size, L~^/^ . This convergence to the 
critical density is indeed visible in the insets of Figs. [Sj a) 
andEJa), although we have not attempted to verify the 
anticipated scaling behaviour explicitly. 

With regard to the behavior of the coexistence num- 
ber densities as a function of parent density ny'i away 
from criticality (Fig. |3J|, we find that the liquid phase 
density varies non-monotonically, while the vapor den- 
sity decreases as the gas phase branch of the cloud curve 
is approached. These trends are also seen in the cor- 
responding MFE results, though the non-monotonicity 
in the liquid density is much weaker and only just per- 
ceptible on the scale of the figure. In the volume frac- 
tion representation (Fig. O, on the other hand, the liq- 
uid volume fraction increases monotonically as the gas 
cloud curve is approached, while the vapor volume frac- 



FIG. 5: (a) MC results for the number densities of the coex- 
isting vapor and liquid phases as a function of n'"^ at T = Tc, 
for (Jc = 1.6. The inset expands the near-critical region and 
shows data for a number of system sizes (expressed in units 
of (j), as indicated, (b) MFE predictions. 



tion decreases. Thus the effect of isothermally reducing 
77,(0) fj-Qj-Q j^g critical value is analogous to that seen on 
lowering the temperature isochorically in a monodisperse 
system: the difference in the volume fraction of the two 
coexisting phases increases. 

This aspect of the phase behaviour is also manifest in 
the "ny' dependence of the surface tension, which may be 
estimated from the form of p{n) via |46j| 



lim ; — T 



In 



pmax^^) 



pmin(„) 



(8) 



where p'^'^'^ {n) / p™™ {n) is the ratio of the peak to trough 
probabilities of the order parameter distribution, which 
provides a measure of the free energy cost of a planar in- 
terface. Formally Eq. ((HJ is valid only in the limit of suf- 
ficiently large system sizes, where the pair of liquid- vapor 
interfaces (whose presence is mandated by the periodic 
boundary conditions) are effectively non-interacting. Al- 
though we have been unable to study systems sufficiently 
large to fully approach this limit in the vicinity of the crit- 
ical point, we do regard our results for 7 as representative 
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FIG. 6: (a) MC results for the volume fraction of the coex- 
isting vapor and liquid phases as a function of n^^' at T = Tc, 
for (Jc — 1.6. The inset expands the near-critical region and 
shows data for a number of system sizes, (b) MFE predic- 
tions. 



of the bulk for parent densities n*^°^ < 0.29. Fig. [7|shows 
our estimates of the surface tension for three isotherms 
having temperatures T < Tc, T = Tc and T > Tc- One 
sees that in each instance the surface tension starts at a 
large finite value on the gas branch of the cloud curve, 
but rapidly decreases with increasing parent density. On 
approaching the high density cloud point it falls to zero 
for T = Tc as expected for critical phase coexistence, 
while for T > Tc and T < Tc it decreases to a finite value 
as is more clearly visible in the inset of Fig. [7| 

Finally in this section, we point out qualitative differ- 
ences in the nature of the coexistence behaviour above 
and below the critical temperature. Fig.|Hlplots our sim- 
ulation estimates of ^, the fractional volume of the liquid 
[cf. Eq. (QJ], as a function of n^^^ for the same three 
isotherms as considered above. For the case T < Tc, the 
fractional volume of the liquid is zero at the gas phase 
cloud, and steadily increases with n^°^ - in a non-linear 
fashion - to reach ^ = 1 at the liquid phase cloud point. 
By contrast for T = Tc, £, increases from zero at the gas 
phase cloud, but reaches a maximum of ^ = 0.5 at the 
critical point. Finally for T > Tc, £, starts from zero, in- 
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FIG. 7: Surface tension 7 for three isotherms at T > Tc, T — 
Tc and T < Tc and size cutoff CTc — 1.6; the inset shows the 
same data on a logarithmic scale to emphasize the behaviour 
near the respective liquid (high density) cloud points. 



creases to a maximum, and then decreases to zero again 
as the high density branch of the cloud curve is reached. 
Thus the distinguishing feature of the "supercritical" (in 
temperature, i.e. T > Tc) coexistence region is that the 
fractional volume of the liquid phase never reaches unity. 
The corresponding MFE predictions show the same qual- 
itative behaviour (data not shown). 
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FIG. 8: MC results for the fractional volume of the liquid 
phase ^ as a function of parent density n*^"^ across the coexis- 
tence region, for three isotherms T > Tc, T = Tc and T < Tc 
and size cutoff CTc = 1.6. Inset; Same data on a logarithmic 
scale. 



IV. PHASE BEHAVIOUR: CUTOFF EFFECTS 

We return now to address the striking feature evident 
from Figs. 121 and 21 of a large difference in the measured 
cloud and shadow curves for the two cutoffs ac — lA and 
1.6. Specifically, the gas phase cloud curve for ac — 1-6 



8 



is shifted to much lower densities compared to that for 
CTc = 1.4. This occurs even though both values of a^ 
are far in the tail of the parent distribution [cf. Fig. [2]. 
The origin of this effect is to be found in the character 
of the particle size distributions in the liquid daughter 
phase. Fig. El shows the form of this distribution for se- 
lected values of n^^^ on the critical isotherm for Oc — 1.6. 
One observes that as the gas cloud point density is ap- 
proached from above, there is a progressive accumulation 
of weight in the large-cr region of the distribution. Thus 
despite the fact that particle sizes around a^. are very 
rare in the parent, they occur (by virtue of fractiona- 
tion) in much higher concentrations in the liquid. The 
physical basis for this is the stronger interaction between 
the larger particles [cf. Eq. 0]; an enhancement in the 
concentration of such particles then yields a substantial 
free energy gain at the shorter interparticle separations 
of the liquid. 
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FIG. 9: Liquid phase size distributions /l(o") for a selection 
of values of n*"' spanning the coexistence region at T = Tc. 
Data are shown for cutoff CTc = 1.6. (a) MC results, (b) MFE 
predictions; densities were chosen to be at the same relative 
positions across the coexistence region as for (a). 

Clearly, therefore, the choice of cutoff has a profound 
effect on the liquid daughter phase distribution, par- 
ticularly in the low density region close to the cloud 
point where the fractionation-induced enhancement of 
the large-cr tail of the size distribution in the liquid 



is greatest. The magnitude of the truncation effect at 
the gas phase cloud point, for T = Tc, is quantified in 
Fig. [TUT a'). which compares the measured forms of the 
shadow phase size distribution at cutoffs a^ = lA,ac = 
1.6 and (Tc — 1.8. For these cutoff values, we find that 
the density at a — ac is enhanced compared to the corre- 
sponding parent density by factors of 11.9(1), 175(7) and 
5280(50) respectively. Also shown, in Fig. Iiur b'l. are the 
corresponding MFE results for the same choice of cutoffs, 
together with a further distribution for the larger cutoff 
tTr = 3; we discuss the form of the latter below. 
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FIG. 10; (a) Size distributions /l(o") in the liquid shadow 
phase at T = Tc for CTc = 1.4, 1.6, 1.8, together with the par- 
ent distribution / (f) plotted for comparison, (b) MFE 
predictions, including the larger cutoff ac — 3. 

Given that the liquid shadow phase distribution is 
highly sensitive to the cutoff and that this phase coex- 
ists with the gas cloud phase, the origin of the sensitivity 
of the locus of the cloud curve to the choice of cutoff [as 
seen in Figs. |31 and 0] becomes clear. We can now also 
rationalize the observation that such significant cutoff ef- 
fects are restricted to parent densities below the critical 
density. For higher densities, the shadow phase is a gas 
of lower density than the parent. In this, the concen- 
tration of large particles is suppressed and that of small 
particles negligibly enhanced because of their weak in- 
teractions. The shadow size distributions are therefore 



concentrated well within the range 0.5 . . .g^ so that no 
noticeable cutoff dependence arises. 

The observed decrease in the gas cloud point den- 
sity n^j with increasing a^ prompts the question as to 
whether the gas phase cloud point density would even- 
tually tend to a zero or nonzero limit as CTc is increased. 
Fig.llllshows the simulation results and MFE predictions 
on the critical isotherm. The former exhibit a further 
strong decrease of n^j by « 25% between Gc = 1.6 and 
1.8; the latter suggest that this trend continues and that 
the cloud point density tends to zero for large Oc- Such 
an unusual effect has previously been seen in theoreti- 
cal studies of polydisperse hard rods with broad length 
distributions 26] and in a Flory-Huggins model for poly- 
mers [43 ■ It is also predicted to occur in solid-solid phase 
separation of polydisperse hard spheres [ig , though only 
for large Uc and distributions with fatter than exponen- 
tial tails. In the present model, however, the decrease of 
nj;[ is clear even for Gc of 0(1), i.e. of the same order 

as a. The physical origin of the decrease of nj;j to zero 
is the appearance (for large CTc) of a second peak in the 
shadow phase size distribution near u^ (Fig. Iior b)). As 
with hard rods, we expect this second peak to eventually 
dominate as Uc increases so that the shadow phase com- 
prises ever more strongly interacting particles whose sizes 
are concentrated near Gc- The appearance of the second 
peak in the size distribution also correlates directly with 
a significant extension of the coexistence region towards 
smaller parent densities, as shown in Fig. ^1 This ex- 
tension shifts the cloud point to lower densities, leaving 
in its wake a region where the fractional volume occu- 
pied by the liquid phase is extremely small, below 10~^. 
This region then crosses over at parent densities around 
•nS^i = 0.06 to more conventional phase behaviour, where 
^ eventually becomes cutoff- independent. These qualita- 
tive features are similar to those observed to hard rods 
with broad length distributions [26|. 
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FIG. 11: The variation of the gas cloud point density njij at 
T = Tc as a function of CTc, as obtained from MC simulations 
(squares) and MFE calculations (circles). 
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FIG. 12: Plot of the fractional volume of the liquid phase 
versus parent density as predicted by the MFE theory at the 
critical temperature t = tc = 0.437. Results are displayed for 
a selection of size cutoffs as shown. 



It is natural to enquire as to the conditions under which 
cutoff dependences can be expected be occur. The an- 
swer will, in general, depend on both the choice of the 
form of the parent size distribution /^'^■'(ct) and the size- 
dependence of the strength of the interaction between 
particles. It is straightforward to show 5] that the den- 
sity distribution of the liquid shadow phase is related to 
the parent via: 



p(")(a)exp[-/3M^!j(a)] 



(10) 



when the cloud point density is small. Thus when the ex- 
cess activity exp[/3^^^(a)] in the shadow phase decreases 
faster with a than the corresponding decay in the parent 
size distribution, the shadow density distribution will be 
peaked at the cutoff. 

We now analyse this scenario in more detail, for a 
more general size-dependence of the interaction strength, 
Cij — {fJiUj)"'. Our hypothesis is that, for suitable par- 
ent distributions /'"'(cr), the cloud point will move to 
vanishing density as the cutoff Cc is made large, while 
the corresponding shadow phase will consist only of par- 
ticles with sizes close to Uc. The shadow is then a 
scaled version of a "standard" (monodisperse, unit par- 
ticle diameter) LJ system, but effectively at a temper- 
ature T/(T^" because all interactions are stronger by a 
factor (crcr')" ~ ^c"! ^^^^ effective temperature decreases 
to zero as Uc is made large. To coexist with the (in- 
finitely, in the limit Cc — > 0) dilute cloud phase, the 
shadow also has to be at zero pressure. At these zero 
temperature, zero pressure conditions, the shadow phase 
will be a crystalline solid, more precisely the ground state 
of the effectively monodisperse LJ system. To check 
that this scenario is self-consistent, we now need to work 
out the excess chemical potential if^{(j) of such a solid. 
It is convenient to use the Widom insertion expression 
exp[— /3/i^x(^)] = (6xp(— /3uo-)); the average is over the 



10 



position where the test particle of diameter <j is inserted, 
uniformly over the whole system, and Uo- is the interac- 
tion potential of this particle with all others. Given the 
assumed size-dependence of the interaction strength, we 
have Ua — (crtTc)"wi, where u\ is the interaction potential 
of a unit diameter test particle inserted into a standard 
LJ solid. Thus exp[-/3/i^h(-^)] ^ (exp[-/3(fTcrc)"ui]). 
When (Jc is large, the average will be dominated by the 
minimum value ui.min < of ui, i.e. the optimal insertion 
positions. Ignoring subexponential corrections then gives 
Mcx('^) = —(fo'c)" 1^1, mini- The cxcess chemical potential 
Q.t a = Oc thus scales as ct^"; looking at Eq. (|10|) . the 
hypothesized cutoff-dominated shadow phase can then 
exist for parent size distributions decaying more slowly 
than /'■^•'(cr) '^ exp(— const x cr^"). Toward smaller a^ 
—fif^{a) ~ (ctCTc)^" decreases by an amount of 0(1) al- 
ready for CTc — CT ^ fc^^" ^ <^c- This fast decrease 
shows that also the assumption of a shadow size distri- 
bution sharply peaked towards CTc is self-consistent. For 
the particular case a = 1 studied in the rest of this paper, 
we have that the limiting parent size distribution where 
cutoff effects start to appear is Gaussian. The Schulz 
distribution has a substantially fatter (exponential) tail 
and it is therefore clear that cutoff dependences should 
arise, as observed. 

As we have seen, the MFE results (cf. Fig. llO() predict 
that, as fJc is increased, so the shadow phase comprises 
ever more strongly interacting particles whose sizes are 
concentrated near <Tc. For large enough cutoffs, the gen- 
eral arguments advanced above show that this shadow 
phase must be a solid because it is effectively at very low 
temperature (as well as low pressure). Since the simple 
approximate free energy used for our theoretical predic- 
tions does not include a branch that could describe such 
solid phases, we have attempted to investigate this possi- 
bility via simulation. Unfortunately, owing to the compu- 
tational burden of simulating systems with large cutoffs, 
we could not access the coexistence regime directly. We 
were, however, able to study the regime of metastable 
coexistence lying between the cloud curve and a system- 
size dependent "effective spinodal" [43. For large cutoffs, 
this region occurs at very low n^°\ and accordingly one 
can, to a good approximation, assign the requisite chem- 
ical potential distribution on the basis of the value of the 
second virial coefficient, as described in Appendix IXI 

Using this approach, we have studied the form of the 
size distribution in the metastable liquid daughter phase 
for a range of large cutoff values. The results in Fig. [TSl 
show that, in accord with the MFE predictions, as (Jc 
increases there is an accumulation of weight at a — (Jc- 
Interestingly too, we find that for CTc — 2.8 the liquid 
spontaneously freezes to an f.c.c. crystal structure, as 
shown in Fig. ^^ We emphasize that this occurs for 
7i(°) values close to the effective spinodal point, not at 
the cloud point itself. However, given the above theo- 
retical considerations it seems likely that, for ac values 
larger than those presently accessible to direct simulation 
at the cloud point, the freezing will occur from the stable 
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FIG. 13: Size distribution /l(o") in the metastable liquid 
phase for cutoffs in the range CTc = 2.0 to CTc = 2.8. Also 
shown for comparison is the parent form /'"'(ct). 



liquid phase 




FIG. 14: Snapshot configuration of the quasi-monodisperse 
metastable solid which coexists with the gas phase at n'"' — 
1 X 10"'^ for ac = 2.8. 

In view of the likelihood of gas-solid coexistence at 
large CTc and small n'-°\ one can speculate as to the 
character of the overall phase diagram in this regime. 
One possibility is that depicted schematically in Fig. [T51 
Here on increasing n^^' from the stable gas region, at 
low T, the system reaches a gas-solid (GS) cloud point 
at which the gas splits off an infinitesimally small quan- 
tity of quasi-monodisperse solid. However, as more of 
this solid is formed it must become more polydisperse in 
order for the system overall to preserve the parent size 
distribution. Hence on further increasing n'"' , the in- 
creasing polydispersity of the solid forces it to split off a 
liquid phase and the system must enter a region of gas- 
liquid-solid (GLS) coexistence. Finally, once all the solid 
has melted, there is a regime of gas-liquid (GL) coexis- 
tence, as is also seen for smaller cutoffs. It seems likely 
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that at higher temperature the soHd phase would not be 
stable, leading to a gas-liquid-solid triple point as shown. 



Triple point 




GS GLS 



FIG. 15: Schematic representation of a possible phase dia- 
gram for large size cutoffs a a- The letters indicate the nature 
of the phases (G: gas, L: liquid, S: solid). 



gas and a quasi-monodisperse crystal. 

Finally, with regard to outstanding questions that have 
not been considered in the present work, we would high- 
light the nature of the critical behaviour. Specifically, 
nothing is known concerning the existence and nature of 
singularities in the cloud and shadow curves in the criti- 
cal region, or indeed what constitutes a suitable choice of 
order parameter for characterizing the critical behavior. 
Intriguing experimental results on polydisperse polymers 
|5l| suggest that the critical exponents are Fisher renor- 
malized, though the reasons for this appear unclear at 
present. We hope to investigate these issues in future 
work. 



APPENDIX A: VIRIAL ESTIMATE OF 

CHEMICAL POTENTIAL DISTRIBUTION AT 

LOW PARENT DENSITY 

For a monodisperse system, the chemical potential can 
be written to first order in the density as 



V. CONCLUSIONS 

In summary, we have deployed state-of-the-art MC 
simulation techniques and MFE calculations to study in 
detail the phase behaviour of a model fluid in which poly- 
dispersity affects both the particle sizes and the strength 
of their interactions. The latter aspect in particular is 
primarily responsible for a dramatic separation of cloud 
and shadow curves compared to a previous study of the 
purely size-disperse case [l^l- We have determined the 
cloud and shadow curves for our model, as well as the 
phase behaviour along certain dilution lines which span 
the coexistence region. We find that the locus of the 
cloud curve is acutely sensitive to the choice of the up- 
per cutoff on the parent particle size distribution, even 
when this cutoff lies far in the tail of the distribution. 
Such effects imply that in experiments on polydisperse 
fluids (see e.g. (ij) it may be important to monitor and 
control carefully the tails of the size (or charge, etc) dis- 
tribution. Otherwise undetected differences could lead 
to large sample-to-sample fluctuations in the observed 
phase behaviour. 

The origin of the observed cutoff dependences has been 
traced to extremely pronounced fractionation effects, the 
nature of which we have elucidated in terms of the char- 
acter of the size distribution of the shadow phase. We 
have also provided criteria for determining which com- 
bination of parent size distributions and interaction po- 
tentials can be expected to display cutoff effects. Our 
theoretical considerations, together with additional sim- 
ulation evidence, suggest that, in the limit of very large 
cutoffs, the cloud point density tends to zero and new 
phases appear, such as a region of coexistence between a 



/3Mcx = 2pa^B*{T) 



(Al) 



with B2{T) the reduced second virial coefficient |53 |. 

For our polydisperse LJ potential, Eq. iPJ, this gener- 
alizes to 



/3Mcx(ff)=2/ da' p{a') 
Jo 



BUT/aa') (A2) 



with 



B^ = -27r 



dr r exp 



-4CTCr' 

T 



0"-0' 



- 1 



Numerical evaluation of the double integral in (|A2|) is 
simplified by noting that for the (untruncated) Lennard- 
Jones potential, i?2(^) can be expressed as a series: 



with r(z) the complete Gamma function. We have not 
corrected for the presence of a potential truncation. 
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